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ABSTRACT 

We present further investigations of the Wilkinson Microwave Anisotropy Probe 
(WMAP) data by means of the Minkowski functionals and the scaling index method. 
In order to test for non-Gaussianities (NGs) with respect to scale-dependencies we use 
so-called surrogate maps, in which possible phase correlations of the Fourier phases of 
the original WMAP data and the simulations, respectively, are destroyed by applying 
a shuffling scheme to the maps. A statistical comparison of the original maps with 
the surrogate maps then allows to test for the existence of higher order correlations 
(HOCs) in the original maps, also and especially on well-defined Fourier modes. 

We calculate the a-normalised deviation between the Minkowski functionals of 
original data and 500 surrogates for different hemispheres in the sky and find eclip- 
tic hemispherical asymmetries between the northern and southern ecliptic sky. Using 
Minkowski functionals as an image analysis technique sensitive to HOCs we find devi- 
ations from Gaussianity in the WMAP data with an empirical probability p > 99.8% 
when considering the \ow-£ range with A£ = [2,20]. The analysis technique of the 
scaling indices leads to the same results for this £ interval with a slightly lower devi- 
ation but still at p > 99.8%. Although the underlying foreground reduction methods 
of the maps differ from each other, we find similar results for the WMAP seven-year 
ILC map and the WMAP seven-year (needlet-based) NILC map for deviations from 
Gaussianity in the low-^ range. Our results point once more to a cosmological nature of 
the signal. For a higher £ range with A£ = [120, 300] the results differ between the two 
image analysis techniques and between the two maps which makes an intrinsic nature 
of the signal on this £ range less likely. When we decrease the size of the analysed sky 
regions for the low-£ study, we do not find signatures of NG in the northern ecliptic 
sky. In the south we find individual spots which show deviations from Gaussianity. 

In addition, we investigate non-Gaussian CMB simulations that depend on the 
/NL-parameter of the local type. These simulations with f^£ al = ±100, ±1000] 
cannot account for the detected signatures on the \ow-£ range. 

Key words: Cosmology, Cosmic Microwave Background, Methods, Non- Gaussianity, 
Minkowski, Observations, Statistical methods, WMAP 
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The primary anisotropies of the Cosmic Microwave Back- 
ground (CMB) are caused by the primordial inhomogeneities 
of the Universe. They are assumed to be generated during 
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a period of cosmic inflation and became the seeds for the 
structure of the density distribution in the Universe we ob- 
serve today. The CMB was emitted from the surface of last 
scattering about 380,000 years after the Big Bang and the 
analysis of its AT/T ~ 10 -5 temperature fluctuations help 
unraveling the mysteries of inflation. 

The WMAP high resolution maps of the CMB allow 
detailed tests on the nature of the primordial density per- 
turbations. The analysis of the CMB power spectrum is al- 
lowing for a high-precision determination of the fundamental 
cosmological parameters. However, any cosmological infor- 
mation that is encoded in the phases and correlations among 
them is not contained in the power spectrum and has to be 
extracted from measurements of higher-order correlations 
(HOCs). 

Evidence for the existence or non-existence of non- 
Gaussianity (NG) is necessary to choose the correct class 
of scenarios for the very early Universe. The inflationary 



scena rio ( |Guth||198l] |Linde||1982| |Albrecht k Steinhardt 



1982) was proposed about 30 years ago to solve the flat- 
ness, monopole and horizon problems of standard Big Bang 
cosmology. It endorses homogeneity and isotropy of the Uni- 
verse, provides a mechanism for the generation of primor- 
dial, approximately scale- invariant and approximately Gaus- 
sian density perturbations and its predictions are consistent 
with the observed power spectrum ( Komatsu et al.||2011 ). 

The Gaussian distribution of scalar (density) and ten- 
sor (metric) perturbations is a generic prediction of single- 
field slow-roll inflationary models. Yet, it has been shown 
that single field as well as two-/multi-field models ( |Linde| 
k Mukhanov| |1997[ |Bartolo, Matarrese k 



Bernardeau k Uzan 2002 ; Vernizzi k Wands 2006 ) generate 



a small amount of non-Gaussianity, below our current ex- 
perimental limits though ( Maldacena|2003 Acquaviva et al. 
2003). Models where the primordial density perturbations 
are generated by a curvaton field may predict a high level of 



NG (Moroi k Takahashi 


2001 


Enqvist k Sloth 2002 Lyth, 


Ungarelli k Wands|2003 


. There are many more alternative 



inflationary scenarios which could generate NG at an ob- 



servable level, e.g. Dirac-Born-Infeld inflation (Silverstein k 



Tong||2004| |Alishahiha, Silverstein k To ng 2004) or ghost 
inflation ( Arkani-Hamed et al.||2004 ). See also review arti- 



cles on NGs from inflationary models written by |Bartolo| 
|et al.| ( |2QQ4| ) and |Chen| ( |2010| ). The determination of non- 
Gaussianity would make it possible to distinguish between 
these different inflationary models. However, inflation still 
remains as a paradigm and a determination of NGs would 
also constrain alternatives, e.g. ekpyrotic and cyclic models 
( Buc hbinder, Khoury k Ovrut||2008] |Lehners k Steinhardt | 
2008 ), and help unraveling the nature of the primordial den- 
sity perturbations in general. 

Various methods of statistical analysis of the WMAP 
data claim the detection of non-Gaussianity and different 
anomalies like hemispherical asymmetries, lack of power at 
large angular scales, alignment of multipoles, detection of 
the Cold Spot, etc. (|Park||2004| |Eriksen et al.||2004[ ~ 



2007) |Hansen, Banday k G6rski||2004 



Hansen et 



de Oliveira-Costa et al.||2004| |Vielva et al. 



manith et al. 



2009) 



aD 



2004; 



2005 



2009 
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In these studies, the level of NG is 
assessed by comparing WMAP data with simulated model- 
dependent CMB maps based on cosmological models and/or 
specific assumptions about the nature of non-Gaussianities 



as parametrised with e.g. the scale-independent scalar pa- 
rameters /nl and #nl • Fewer studies are testing the random 
phase hypothesis for Gaussian random fields analysing the 
distribution of the Fourier phases of the data \ Chiang et al. 
20031 |Coles et al.|2004| [Naselsky et al.|2005| [Chmng, Nasel- 
sky k Coles|2007| ). These model-independent tests also show 
signatures of anomalies and non-Gaussianities. 

In this paper, we apply the well established method of 
surrogate data sets to analyse Fourier phase correlations in 
the WMAP data. The method of the surrogates has first 
been established in Theiler et al. (1992) in order to de- 



tect weak nonlinearities in time series. For given possibly 
nonlinear data, so-called surrogate data sets are generated, 
which mimic the linear properties of the original time series. 
A comparison of original data set and surrogates through 
statistical measures sensitive to higher order correlations 
(HOCs) can reveal a significant deviation from linearity in 
the data. Extensions to this formalism to three-dimensional 
galaxy distributions ( Rath et al.|2002" ) and two-dimensional 
simulated flat CMB maps (Rath k Schuecker 2003) have 



been proposed and discussed. In a prior report, we intro- 
duced a two-step surrogatisation scheme for full sky CMB 
observations which allows for a scale-dependent analysis of 
the data (Rath et al. 2009). This method yields significant 
signatures for both NG and ecliptic hemispherical asymme- 



tries (Rath et al. 2009[ 2011), especially on largest scales 



(£ ^ 20). In these papers, the HOCs are measured with the 
weighted scaling index method (SIM) introduced in Rath 
et al.|(|2002|),|Rath k Schuecker | ( pOO^ and |Rath, Schuecker 
k Banday| ( |2007| ). The image data is represented as a point 
distribution comprising the spatial and scalar information 
of each pixel. The scaling indices then offer the possibility 
to estimate the local scaling properties of such a point dis- 
tribution and depend on a scale parameter r. 

We concentrate on the set of statistics known as 
Minkowski Functionals (MFs) as a comparative measure 
sensitive to HOCs. A general theorem of integral geometry 
states that three Minkowski Functionals quantify the inte- 
grated geometry and topology and therefore morphology of 
a two-dimensional density field. For a pixelised temperature 
map AT/T(n), we consider the excursion sets of the map, 
defined as the set of all map pixels with value of AT/T(n) 
greater than some threshold level v. The Minkowski func- 
tionals of these excursion sets completely describe the mor- 
phological properties of the underlying temperature map 
AT /T(n). These measures embody the information from all 
orders of the correlation functions. They are additive mea- 
sures which allows us to calculate them efficiently by sum- 
ming up their local contributions. The calculation of n-point 
correlation functions is much more expensive computation- 
ally ( |Eriksen etld~1|2005| |Szapudi||2009| ). In comparison to 
the SIM the Minkowski functionals are scale-independent, 
i.e. they do not depend on any scale parameter. 

We apply our methods to WMAP seven-year ILC and 
NILC maps. Furthermore, we simulate and analyse non- 
Gaussian CMB maps that depend on the scale-independent 
parameter and serve as first toy models for further 
tests on inflationary models. 

The paper is organised as follows. In Section 2 we briefly 
describe the observational and simulated data we use in our 
study. The method of generating surrogate maps is reviewed 
in Section 3. A discussion and comparison of the used test 
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statistics, Minkowski functionals and scaling indices, can be 
found in Section 4. In Section 5 we present our results and 
we draw our conclusions in Section 6. 



2 DATA SETS AND SIMULATIONS 

For our studies we use the WMAP seven-year foreground- 
cleaned internal linear combination map (in the follow- 



ing: ILC7) generated and provided by Gold et al. (2011). 
This map can be directly do wnloaded from the LAMBDA- 
websitej^ with a HEALPhQ ( |G6rski et aLl|2005| ) resolution 
parameter of A s id e = 512, corresponding to 3, 145, 728 sky 
pixels. The ILC7 map has one degree FWHM. For compar- 
ison we also include the seven-year needlet-based ILC map 
(in the following: NILC7) produced by |Basak Delabrouille 
poll) pursuing a different approach for foreground removal 
Because the needlet implementation of the ILC allows for op- 
timising weights both as a function of sky direction and as 
a function of scale, the NILC7 is significantly less contam- 
inated by foregrounds than other existing maps obtained 
from WMAP, in particular at low Galactic latitude (see|De-| 



labrouille et al. ( 2009 ) for a description of the needlet ILC 



method). It features a better total resolution, correspond- 
ing to the WMAP W-band resolution, and has an original 
HEALPix resolution of A side = 1024. We downgrade the 
map to A s ide = 512 in the employed HEALPix software. The 
ILC7 as well as the NILC7 map are weighted linear combina- 
tions of the five frequency channels K, Ka, Q, V and W that 
recover the CMB signal even in heavily foreground affected 
regions. The weights are calculated by requiring minimum 
variance in a given region of the sky under the constraint 
that the sum of the weights is unity. Full sky maps ensure 
the required orthogonality of the set of basis functions Ytm 
when generating surrogates. 

For our comparison of the ILC7 and NILC7 maps on 
larger scales, 120 ^ £ ^ 300, we choose to decrease the reso- 
lution of the original NILC7 map down to one degree FWHM 
before downgrading to a HEALPix resolution of A s id e = 512. 
In these higher ^-ranges the beam effects have stronger in- 
fluences. For our studies of the low-^ ranges with ^ £ ^ 20 
this influence is negligible. 

The primordial non-Gaussianity that might arise dur- 
ing cosmic inflation has often been parametrised by the /nl 
parameter in the following simple form of the curvature per- 
turbation $ (Bardeen|1980) with quadratic correction 



$(x) = ®l{x) + /nl^nl(x) 



(1) 



where <i>z, are Gaussian linear perturbations and $nl(x) is 
defined as 



(2) 



and /nl characterises the dimensionless amplitude of 
quadratic correction. This form of /nl describes non- 
Gaussianity of the local type. In order to test for possible 
impacts of local type NG on the data, we compare the re- 
sults obtained from the original data to models based on a 
scale-independent /nl^ parameter. Therefore, we compute 



1 http://lambda.gsfc.nasa.gov 

2 http://healpix.jpl.nasa.gov 



temperature maps as constrained realisations of non-linear 
fields with varying /nl 3,1 . 

The spherical harmonic coefficients aim of the CMB 
temperature (i — T) and polarisation (i — E) anisotropics 
in harmonic space are related to the primordial fluctuations 
®im(x) via 



/ 



diXX 2 ®£m(x) Ot\(x) 



(3) 



as a function of comoving distance x, where a\(x) is the real 
space transfer function. A set of linear and non-linear spher- 
ical harmonic coefficients are a valid realisation of temper- 
ature and polarisation fluctuations for a given cosmological 
model, in this case ACDM model. The total aim coefficients 
and a desired level of f^ al NG are then calculated as 



linear . p non — linear 

Mm + /nl * a em , 



(4) 



where the linear and non-linear a£ m coefficients are provided 
by |Elsner fe Wandelt| fl2009| ). 

To explore the plausibility of /nl as an explanation for 
the observed large-scale NGs, we use these a^ m coefficients 
to simulate five co-added VW-band maps with f^£ al = 



[0, ±100, ±1000]. As shown in |Rath et aL] ( |2Qll| ) the results 
for the method of surrogates of the simulated VW-band 
maps are similar to those of ILC-like simulations on large 
scales. The specific beam properties of the WMAP satellite 
are taken into account using the respective window func- 
tions for each differencing assembly (VI — V2, Wl — W4), 
being available again on the LAMBDA- website 1 . For every 
assembly, we add Gaussian noise to these maps with a par- 
ticular variance for every pixel of the sphere. This variance 
depends on the number of observations Ni(0,(f)) in the re- 
spective direction and the noise dispersion per observation 
per different assembly, gq^. After this procedure, we accu- 
mulate the V- and W-bands to a co-added VW-band via a 
noise- weighted sum ( |Bennett et al. 2003): 



T{9,< 



(5) 



In this equation, A characterises the set of re- 
quired assemblies, for the co- added VW-map A = 
{VI, V2, Wl, W2, W3, TV4}. The parameters = [0, tt] and 
cj) = [0, 2n] correspond to the co- latitude and the longitude 
on the sphere, while the seven-year noise per observation of 
the different assemblies is denoted by <jo, given by Jarosik 
et al.| ( |201l| ). Again, 
resolution of A s id e 



the maps are decreased to a HEALPix 
= 512. Finally, we remove the residual 



monopole and dipole by means of the appropriate HEALPix 
routine. 



3 METHOD OF SURROGATES 

In order to constrain possible phase-correlations of the 
Fourier phases of the cosmic microwave background, we 
generate so-called surrogate maps as proposed in |Rath fc] 
|Schue cker (2003). We destroy possible phase correlations in 
the data by applying a shuffling scheme to the phases. In or- 
der to test for the scale-dependence of non-Gaussianity this 
is done with a two-step procedure. 

The full sky CMB temperature anisotropy map of the 
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celestial sphere with AT/T(n) = AT/T(9, <fi) can be ex- 
panded in orthonormal spherical harmonics Yg m as 

oo I 

AT/T(9,4>) = J2 Yl "tmY tm (0,<l>) (6) 

1=0 m=—l 

with the complex spherical harmonic coefficients 



j dn T(n)F/ m (r 



(7) 



where n is the unit direction vector and aim = \agrn\e 1 m • 
The linear properties of the underlying density field are con- 
tained in the absolute values |a^ m |, whereas all HOCs - if ex- 
istent - are encoded in the phases <j)£m and the correlations 
among them. If the spherical harmonic coefficients are inde- 
pendent Gaussian random variables their probability density 



P(atm)dain 



da£ n 



(8) 



(e.g. Rubakov & Vlasov (2010J) and the ortho normality re- 
lation holds: 



(9) 



where Ci is the angular power spectrum and 5 is 
the Kronecker delta. In this case, the amplitude \at m \ 
would be Rayleigh-distributed and the phase <p£m — 
arctan (Im(ag m )/Re(ag m )) would be independent and iden- 
tically distributed (i.i.d.) and follow a uniform distribution 
in the interval [— 7r,7r]. The temperature values themselves 
would be normally distributed as well. The measured angu- 
lar power spectrum C° bs depends on the a£ m coefficients by 



Cobs 
P — 



2£ + l 



(10) 



where (C^ bs ) = Cg. In case the CMB variables ae m are in- 
dependent and Gaussian distributed and their phases are 
therefore i.i.d. (independent and identically distributed) and 
consequentially uncorrelated, it is adequate (iKomatsu et al. 



(2009) and references therein) to only investigate the lin- 
ear properties of the data described by the power spectrum 
Cg. In this paper, we focus on the non-linear information 
and test for possible phase correlations of <j)£ m in the data, 
which we define as non-Gaussianities of the CMB, to check 
if the above requirements are fulfilled. 

To test the hypothesis of independent Fourier phases 
we try to exclude further features of the data as e.g. arte- 
facts due to experimental constraints. Therefore, we make 
sure that the data meets the following two requirements be- 
fore generating the surrogate maps. The temperature distri- 
bution is Gaussian and the set of phases is uniformly dis- 
tributed in the interval [— 7r,7r]. To fulfil these conditions 
we perform two preprocessing steps. First, the almost Gaus- 
sian distributed temperature values of the original map in 
real space are replaced by an exact Gaussian distribution 
in a rank-ordered way, i.e. the lowest value of the original 
distribution is replaced with the lowest value of the Gaus- 
sian distribution etc. Second, in order to guarantee that the 
phases 0£ m are identically distributed the Fourier phases 
are remapped on to a set of uniformly distributed ones. No 
significant dependence on the specific Gaussian or uniform 



realisation, respectively, was found in these preprocessing 
steps. 

The scale dependent surrogate maps that are later 
analysed are obtained as follows. In our scale-dependent 
analysis we focus on the two ^-ranges A£i = [2,20] and 
A£ 2 — [120,300]. The interval A£i covers the largest spa- 
tial scales of the CMB, while the first peak of the power 
spectrum lies in the interval A£ 2 . Our previous results have 
shown that findings of an almost scale-independent analysis 
with a shuffling range Ai = [2, 1024] are only the superposi- 
tion of signals from the analysis with A£± and A£ 2 shuffling 
ranges ( |Rath et al.||201l| ). Therefore, we assume that the 
results of these two intervals are of special interest in the 
analysis of phase correlations. We first generate a first order 
surrogate map, in which any correlation of phases <^ m with 
£ outside the certain A^-range of interest is destroyed: The 
phases <^ m with 2 < £ < 1024 and £ g A£ = [£min,£max], 
< m < £, are randomised through a shuffling procedure. 
In a second step, N (N = 500 throughout this study) real- 
isations of second order surrogate maps are generated from 
the first order surrogate map: The remaining phases <^ m 
with £ G A£ = [£ m in,£max] are shuffled, while the already 
randomised phases in the first order surrogate for the scales 
not under consideration are preserved. Note that the Gaus- 
sian properties of the maps, which are given by |a^ m |, are 
exactly preserved for all surrogate maps. 

The first order surrogate preserves phase correlations 
in the A£ range, if they were present in the original maps. 
In the second order surrogate map, where all phases are 
randomly distributed, all original correlations will be com- 
pletely destroyed. The statistical comparison of the two 
classes of surrogates will thus reveal possible higher order 
correlations in the original maps on defined Fourier modes 
and uncover signatures for deviations from Gaussianity. 

One might argue that the above discussed phase shuf- 
fling is not the intuitive approach, since it is also possible to 
construct surrogates (of first as well as of second order) by 
just replacing the respective Fourier phases with randomly 
generated ones. In order to account for this aspect, we addi- 
tionally investigate surrogates (for selected data sets) that 
were constructed by replacing the phases with a completely 
new set of uniform distributed values in the interval [— 7r, tt]. 
Note that for this replacement approach, one has to generate 
a new set of phases for the relative ^-interval. It might be 
preferable to only rely on the information which is given in 
the underlying data set, instead of adding something new. 

Before analysing the maps with different image analysis 
techniques we decrease the HEALPix resolution of the maps 
from iV s ide = 512 to N s id e = 256 for the data analysis and 
^Vside = 64 for the simulated maps. 



4 TEST STATISTICS 

4.1 Minkowski Functionals 

A full morphological specification of an image requires ge- 
ometrical as well as topological descriptors to characterise 
not only the shape and content but also the connectivity of 
spatial patterns. Hadwiger's theorem (Hadwiger 1957) in the 
field of integral-geometry states that any complete morpho- 
logical descriptor of a set Q in d-dimensional Euclidian space 
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E is a linear combination of only d + 1 functionals which 
meet some simple requirements; the so-called Minkowski 
Functionals Mj with j ranging from to d. Functional Mo 
and Mi in two-dimensional space can be interpreted as the 
familiar geometric quantities surface area and perimeter, re- 
spectively. Functional M2 is the topological Euler charac- 
teristic. In spherical, two-dimensional CMB sky map space 
S 2 of radius R and a constant curvature K = R~ 2 they are 
formally defined as 



ing an edge or vertex more than once. As suggested by 



Mo 



Mi 



M 2 



f da 

Jq 

4JdQ 

2^ L 



dl 



di h Q 



where da and dl denote the surface element of S 2 and the line 
element along the smooth boundary dQ, respectively, as ex- 
plained in |Schmalzing &; Gorski| p998| . The factor k g is the 
geodesic curvature. The Minkowski functionals have been 
introduced into cosmology as descriptors for the morpho- 



logical properties of large-scale structure by Mecke, Buchert 
fe Wagner I (|1994|) and of CMB sky maps by jWinitzki fe 



Kosowsky (1998) and Schmalzing fe Gorski (1998) 



In order to study the morphology of the temperature 
anisotropics AT/T(n), which can be considered as a smooth 
scalar field on S 2 , we calculate the three Minkowski Func- 
tionals (MFs) of excursion sets Q v in the pixelised spherical 
two-dimensional maps of original data, surrogates and sim- 
ulations. Q u is defined as the set of all map pixels with value 
of AT /T greater than or equal to some temperature thresh- 
old v by Q v = {n e § 2 |AT/T(n) ^ v}. 

The maps we are analysing are pixelised according 
to the HEALPix pixelisation scheme. HEALPix produces 
a partition of a spherical surface into exactly equal area 
quadrilateral pixels of varying shape which simplifies the 
calculation of the Minkowski functionals. The pixel size de- 
pends on the HEALPix resolution parameter of the grid 
equal to iV s id e = 1,2,4,8, ... corresponding to a total num- 
ber of pixels of Afpi x = 12 x iV s 2 ide = 12,48,192,768,... In 
this work we use a resolution parameter of iV s id e = 256 and 
A^ S ide = 64, respectively. 

The temperature maps we want to analyse are divided 
into an active and a non-active part by running over 200 
threshold steps Vi with — 4<jt ^ v% ^ +4<tt- At the first 
threshold step z/o nearly every pixel is included in the active 
part besides few outliers with T < —4<jt- The last step ^199 
excludes most of the pixels. 

We adapted an algorithm of Michielsen &; De R aedt 



( 2001 ) to compute the Minkowski functionals of the pixelised 



maps. Conceptually, each active pixel is decomposed into 4 
vertices, 4 edges and the interior of the pixel. We count the 
total number of active squares n s and edges n e and vertices 
n v between active and non-active pixels and compute the 
area Mo, the integral mean curvature or perimeter Mi and 
the Euler characteristic M2 from 

M n s 

Mi = -4n s + 2n e 
M 2 = n s - n e + n v . 

A technical difficulty with this procedure is to avoid count- 



Michielsen &; De Raedt (2001) we build up the original im- 



age by adding active pixels to an initially empty temporary 
image one by one. Depending on whether the surrounding 
pixels have already been activated in the temporary image 
or not, we then add up edges and vertices to their total 
numbers. The number of arithmetic operations required to 
compute M , Mi, and M 2 scales linearly with the number 
of active pixels and the total number of pixels of the image. 

4.2 Weighted Scaling Indices 

For comparison we assess possible scale-dependent NGs in 
the CMB with the scaling indices method (SIM) as done 
in our series of earlier papers (Rath et al.|[2002 Rath & 



Schuecker 2003; Rath, Schuecker & Banday 2007; Rath et al. 
2009 ;|Rossmanith et al.|2009| IRath et al.|2011| IRossmanith 



et al.|20lT ). The basic idea for this test statistic comes from 
the calculation of the dimensions of attractors in non-linear 
time series analysis |Grassberg er &; Procaccia]|1983| ). The 
scaling indices have been extended to the field of image pro- 
cessing for texture discrimination. If the image data is repre- 
sented as a point distribution in a d- dimensional embedding 
space the scaling indices represent one way to estimate the 
local morphological properties of this point set. Point-like, 
ring-like and sheet-like structures can be discriminated from 
each other and from a random background. 

In order to apply the SIM on the spherical CMB 
data, we have to transform the temperature anisotropics 
AT/T(m) with its pixels at positions m = (0*, i = 
1, Npi X on the sphere to a three-dimensional point distri- 
bution of Npi X points P = {p-} in an artificial embedding 
space, for which then the local scaling properties are assessed 
with the SIM. 

For each point the local weighted cumulative point dis- 
tribution p is calculated as 



(11) 



where s r («) denotes a shaping function depending on a scale 
parameter r and a distance measure In principle any 

differentiate shaping function and any measure can be used 
for calculating the scaling indices. We use the Euclidian 
norm as distance measure and set of Gaussian shaping func- 
tions, which leads to 

Npix _d j± ^ 

p(Pn r ) = e ~~ = llPi-PjII- (l 2 ) 

The weighted scaling indices a(p i ,r) are then obtained by 
calculating the logarithmic derivative of p(p i ,r) with respect 
to r: 



d\og(r) 



(13) 



With the definition in Equation [12] the weighted scaling in- 
dices are expressed by 



(14) 



3 = 1 e 
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4.3 Statistical Interpretation 

The two image analysis techniques, Minkowski functionals 
and scaling indices, are applied to scale-dependent full sky 
surrogate maps and 768 overlapping hemispherical maps, 
with two different solid angles: a solid angle of 2ty (apex 
angle tt) and ~ 0.67T (apex angle tt/2). These hemispheri- 
cal maps are rotated around the full sky to study possible 
scale-dependent phase-correlations in the Fourier space of 
the maps with a certain spatial localisation of the phenom- 
ena. In this paper we focus on the analysis of the Minkowski 
functionals and compare our results to the scaling index 
method. 

In order to quantify the degree of agreement between 
the surrogates of different orders with respect to higher or- 
der correlations found with the three Minkowski functionals 
Mo , Mi and M2 we calculate the mean of each Minkowski 



functional M* 



ro2 for N — 500 realisations of the second 



order surrogate per hemisphere h and threshold bin v 
1 N 

A := (M* jSnrro2 (z/,/l)) = — M ^,surro2(m,U,h) 



N 

and the standard deviation 



0~ M+ 



2 (zy, h) 



( 1 N 

\ m — l 



2(ra, z/, h) — a y 



1/2 



for m = 1, N. We combine mean and standard deviation 
in a diagonal x 2 statistic per hemisphere h for the surrogates 
of first order 



2 



3=0 



L (i/, h) 



,(zy, h) 



and for the second order surrogate maps 
M* ,surro2 (m, v, h) 



2 

Xm. 



2 (h,m) 



3=0 



3 (l/, /l) 



Finally, the degree of agreement between the two types of 
surrogates is quantified by the cr-normalised deviation S 



S(xmA H )) 



(15) 



, (/i)) and 



for each of the 768 hemispheres with (xl 
a^2 (h) denoting the mean and the standard devia- 

XM *,surro2 

tion of xif„ surro2 W f° r the ^ = 500 second order surro- 
gates. Thus, we obtain the a normalised hemispherical devi- 
ations S(x 2 M (h)), S(x 2 Ml (h)) and S(xm 2 W) for tne Ar ea, 
Perimeter and Euler characteristics. 

As for the Minkowski functionals, we calculate the mean 
(a(r)) and standard deviation a a ( r ) of the scaling indices 
a(p-,r) for the set of 768 hemispherical maps. The scaling 
indices are calculated with a scale parameter r = 0.25, which 
corresponds to no as defined in our earlier works, e.g in 
Rossmanith et al. ( 2009 ) . The differences of the two classes 
of surrogates are again quantified by the cr-normalised devi- 
ation S: 

g(Y^ Ysurrol (Xs-urro2) (16) 



where Y represents a diagonal x 2 statistic 



E 

i=i 



Bj - (Bj) 



(17) 



as a combination of the mean and the standard deviation 
where B\(r) — (a(r)), B 2 (r) — cr a (r)- See|Rath, Schuecker 



|&; Banday] ( 2007 ) for a detailed description regarding the 
scaling indices statistics. 

The results we obtain with the hemispherical study of 
the sky are visualised in so-called S'-maps. Each pixel centre 
of a full sky map with a HEALPix resolution of Af s ide = 8 
marks one of the 768 hemispheres' poles. After calculating 
the deviation S between one first order surrogate and 500 
realisations of second order surrogates for each individual 
hemisphere, we plot the respective value in a sky map at 
that pixel position where the z-axis of the rotated hemi- 
sphere pierces the sky. This is done for both, the Minkowski 
functionals and the scaling indices. For one single S'-map 
with 768 values stemming from the comparison of 500 sur- 
rogates of second order with one first order surrogate we 
need to calculate 384,768 hemispherical maps. 



5 RESULTS AND DISCUSSION 
5.1 ILC and NILC Maps 

Figure [l] shows the S- values per hemisphere for a deviation 
from Gaussianity in ILC7 and NILC7 data on largest scales 
with A£ = [2,20], found with the method of surrogates by 
using a shuffling approach and a phase replacement proce- 
dure, respectively, and analysed by the Minkowski function- 
als as well as the scaling index method. The x 2-s tatistics 
of the two image analysis techniques yield consistent results 
on largest scales and show significant signatures for eclip- 
tic hemispherical asymmetries and non- Gaussianity in CMB 
sky maps. The signal for the Minkowski functionals is maxi- 
mal in the southern ecliptic sky whereas for the SIM we find 
the maximum in the northern ecliptic sky. In Table [l] and [2] 
we summarised the deviations S and the empirical probabil- 
ities p of the three Minkowski functionals and the SIM. The 
results are shown for the full sky and a pair of hemispheres 
which consists of the hemisphere with maximum S(x 2 ) an d 
its opposing hemisphere on the other side of the sky. 

The two CMB maps, ILC7 and NILC7, are different 
enough in their implementation that the presence of residual 
foregrounds and noise, if important in the present analysis, is 
not expected to result in similar NG detections. The needlet 
ILC is, in principle, less contaminated by foregrounds and 



noise, as shown on 5-year data by Delabrouille et al. (2009) 



However, the pattern in the S'-maps of these two maps are 
almost identical to each other. The cr-normalised deviations 
S are higher for the NILC map and range up to 9.97 for the 
perimeter with an empirical probability above 99.8%. The 
stronger signatures for NGs in the NILC map can further 
be attributed to the fact that this map did not undergo the 
one-degree smoothing of the ILC7 map. 

Expectedly, the replacement of the original phase distri- 
bution by a distribution of random values in NILC7 during 
the surrogates generating process yields the same pattern 
in the S'-maps as for the shuffling procedure. However, the 
shuffling approach avoids any dependence on additional data 
and is therefore preferable for generating surrogates. 

We show that the scale-dependence of the scaling index 
method is not a limitation to that measure since the scale- 
independent x 2 statistic of the Minkowski functionals still 
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Figure 1. Deviations S^X 2 ) of Minkowski Functionals Mo, M\ and M2 of the rotated hemispheres for the ILC7 (upper row) and NILC7 
map (middle row) in Galactic ccordiantes. In the lower row we show the results of the phase replacement method for NILC7. The plots 
to the very right show the corresponding results S(x 2 ) for the respective maps gained by the scaling index method. The Grange for the 
method of the surrogates is At = [2, 20]. 





Full Sky 


hemisphere 

Smax 


hemisphere 
Opposite S 




Full Sky 


hemisphere 

Smax 


hemisphere 
Opposite S 


x 2 


(S\%) 


(S\%) 


(S\%) 


x 2 


(S\%) 


(S\%) 


(S\%) 


Area 

Perimeter 
Euler 


0.62 | 86.4 
0.93 88.6 
1.44 92.2 


6.72 | 99.6 
7.33 | >99.8 
7.24 >99.8 


3.05 | 98.8 
4.52 99.4 
3.62 99.0 


Area 

Perimeter 
Euler 


1.03 | 88.2 
0.89 86.4 
0.77 84.4 


9.51 | >99.8 
9.97 >99.8 
9.50 >99.8 


5.98 | 99.8 
7.31 99.8 
7.22 | >99.8 


SIM 


0.41 | 57.0 


8.9 | >99.8 


6.1 | 99.8 


SIM 


0.29 | 51.4 


7.53 | >99.8 


6.23 | >99.8 



Table 1. Deviations S calculated for the ILC7 map at a shuffling 
range of Ai = [2,20] from Equation |15| and the empirical prob- 
abilities p of the scale-independent diagonal x 2_s tatistics for full 
sky functionals and two selected single hemispheres for the three 
Minkowski Functionals and the scaling index method. This table 
corresponds to the first row of Figure ^ 



yields the same results as the scale-dependent x 2 statistic of 
the scaling indices, that depends on scale parameter r. 

All three Minkowski functionals area, perimeter and Eu- 
ler as well as the SIM detect phase correlations and therefore 
non-Gaussianities in the data with almost identical spatial 
signatures. However, the single values in the *S-maps do not 
demonstrate local NGs but must be interpreted as an overall 
signal per hemisphere. Note that in this hemispherical study 
large overlapping sky patches have been analysed. S- values 
in the upper sky also account for signal from the lower sky 
and the other way around. 

In Figure |2] and [3] the three Minkowski functionals area, 
perimeter and Euler characteristic of the ILC7 and NILC7 
map, respectively, are plotted as a function of threshold val- 
ues Vi. We compare one surrogate of first order with 20 re- 
alisations of second order surrogates generated for the i- 
range of interest, here A£ — [2, 20]. The full sky Minkowski 



Table 2. The same as Table^ but for the NILC7 surrogate maps. 
This table corresponds to the second row of Figure^ 



functionals do not show differences between the two classes 
of surrogates. Single hemispheres though show clear devi- 
ations between first and second order surrogate maps. For 
the plots we choose again the pair of pixels as defined above: 
the hemisphere with the highest 5'(x 2 )-value and its oppos- 
ing hemisphere. In general, we refer to hemispheres mainly 
located in the northern Galactic sky as upper hemispheres, 
southern sky hemispheres are called lower hemispheres. 

The first order surrogate functionals of the upper or 
lower hemisphere differ from the respective second order 
surrogates for most of the threshold levels v. These devi- 
ations between the two classes of surrogates lead to the as- 
sumption that phase correlations of Fourier phases (j)i m with 
2 ^ / ^ 20, manifested as certain morphological structures 
in the temperature distribution in real space, are existent 
in the first order surrogates, in which only phases outside 
Ai — [2, 20] are randomised, but are destroyed in the sec- 
ond order surrogate maps. 

The first order surrogates of area, perimeter and Eu- 
ler functionals show contrary behavior comparing the upper 
and lower sky. Where the surrogates of first order of the 
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Figure 2. Minkowski Functionals Mo (left), Mi (middle) and M2 (right) of the first (black) and 20 second order surrogate maps (green), 
calculated for A£ = [2, 20]. The surrogate maps were derived from the ILC7 map. The upper row shows the Minkowski functionals for 
the full sky, the middle row for an ecliptic northern sky hemisphere and the lower row for an ecliptic southern sky hemisphere as selected 
from the S'-maps. The coordinates are given in Table [l] 
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Figure 3. Same as Figure [2] but for the NILC7 surrogate maps. 
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Figure 4. Probability density P(a) of the first and second order 
surrogates for the scaling indices calculated for surrogate maps 
with the t interval At = [2,20]. The black lines denote the first 
order surrogates derived from the NILC7 map. The green lines 
are the densities for the corresponding 20 realisations of second 
order surrogates. 



chosen upper sky hemisphere lie above the second order sur- 
rogates it is the lower sky first order surrogates' Minkowski 
functionals that lie below. The algebraic sign of the devia- 
tions between the two classes of surrogates depends on the 
analysed sky region. The results of NILC7 in Figure [3] show 
in principle the same deviations of the two surrogates and 
between northern and southern sky. The absolute amplitude 
of the Minkowski functional M2 (Euler) for the NILC7 map 
though is larger for negative thresholds compared to positive 
thresholds. In the case of the ILC7 map the two amplitudes 
are nearly equal. 

Figure [4] shows the probability densities P(a) of the 
scaling index method of one first and 20 second order sur- 
rogate maps for the ^-interval A£ = [2,20], again for the 
full sky analysis as well as for two opposing hemispheres 
for comparison. The density distributions of the second or- 
der surrogates with fully randomised phases are shifted to- 
wards lower (higher) a values for the upper (lower) hemi- 
sphere. This shift must be interpreted as a global trend in- 
dicating that the first order surrogate map has less (more) 
structure than the respective set of second order surrogates. 
As these deviations of different algebraic sign cancel each 
other for the full sky, we do no longer see significant dif- 
ferences in the probability densities in the full sky analysis. 
The SIM applied to the NILC7 map shows again that the 
the morphology of the temperature distribution depends on 
the analysed sky region and differs in algebraic sign between 
northern and southern Galactic sky, as was the case with the 
Minkowski functionals. The deviations between first and sec- 




euler 



euler 




Figure 5. Significance S per hemisphere between surrogate of 
first and second order of NILC7 calculated for one threshold value 
Vi. For perimeter, S is calculated at the threshold value where 
Mi a 

urroi becomes maximal (upper left). For Euler, we chose 
the threshold values where M2, S urroi becomes minimal (lower 
left) and maximal (lower right), respectively. For comparison, we 
show the deviations S((a)) for the mean of the scaling indices 
(upper right). The ^-range for the method of the surrogates is 
A£ = [2,20]. 



ond order surrogates demonstrate the existence of HOCs in 
the WMAP Fourier phases. In our earlier works we found 
similar results for the ILC7 and NILC5 (five-year needlet- 
based ILC) map (see Rath et al. ( 2011| ). 

The x 2 results discussed above do not depend on the 
algebraic sign anymore. In order to visualise the sign- 
dependence we plot the deviation S per hemisphere be- 
tween first and second order surrogates calculated with the 
Minkowski functionals for only one chosen threshold, shown 
in Figure [5] and keep the algebraic signs of the deviations 
with this calculation. The chosen thresholds correspond to 
the minima or maxima of perimeter and Euler functional. 
For the SIM we choose to plot S(Y) with Y — (a) corre- 
sponding to Equation [16] to keep the algebraic sign, as was 
done in previous works. For the Euler characteristic we see 
that the sign of the deviation depends on the threshold level 



An analysis of different ^-ranges was also performed in 
|Hansen et al. (|2009|. In their work, the dipole directions of 



the power distributions of different multipole blocks of 100 
multipoles each were calculated. It was shown that for the 6 
multipole blocks in £ = [2, 600] the dipole directions lie very 
close to each other. We compare our non-linear statistics for 
the method of surrogates with A£ = [2, 20] with these linear 
findings in Figure [6] We show our coordinates of the hemi- 
spherical pairs defined above and add the directions of the 
6 dipoles of the power distributions. Note that the linear 
and non-linear statistics can have different scale sensitivity. 
Although it is not clear to what extent the results can be rec- 
onciled in detail, interestingly, our large-scale investigations 
with Minkowski functionals and scaling indices as well as the 



results of Hansen et al. ( 2009 ) show evidence of asymmetry. 



As mentioned above, the original power spectrum is exactly 
preserved for all our surrogate maps. The dipole direction 
of the multipole block £ = 2 — 100 lies close to our large- 
scale SIM results for ILC7 and close to the Minkowski (Area, 
Perimeter) results for NILC7. Furthermore, these three pix- 
els are close to the southern ecliptic pole. 
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perimeter 



Figure 6. A sky map in Galactic coordinates with the directions 
of maximal signal for NGs found in ILC7 (red) and NILC7 (blue) 
with Minkowski functionals (Area, Perimeter: blue, Euler: dark 
blue, Area, Perimeter, Euler: dark red) and SIM (light blue, light 
red). The yellow pixels show the dipole directions of the power 
distribution in 6 blocks of 100 multipoles each taken from [Hansen 
|eTal^ ( |2009] ). The green pixels correspond to the ecliptic northern 
and southern pole, respectively. 




Figure 8. Deviations S(x 2 ) of the Minkowski Functionals Mo, 
Mi and M2 and the scaling indices (from left to right and top 
to bottom) of the rotated hemispheres for the NILC7 map for 
the smaller apex angle tt/2. The ^-range for the method of the 
surrogate is A£ = [2, 20]. 




Figure 7. Deviations of the Minkowski Functionals Mo, Mi and 
M2 of the rotated hemispheres for the ILC7 (left column) and 
NILC7 (right column) map. The ^-range for the method of the 
surrogates is Al = [120,300]. 



In Rossmanith et al. (2012) we introduced the method 



of surrogates for an incomplete sky. We find that even when 
the complete Galactic plane is removed, NGs and hemispher- 
ical asymmetries can still be detected in the CMB and con- 
clude that the Galactic plane cannot be the dominant source 
for the found anomalies. Our results point in general towards 
a violation of statistical isotropy of the Universe and disfavor 
single-field slow roll inflation. 

Supplemental to the results on largest scales, Figure [7] 
shows the deviations S(x 2 ) P er hemisphere of the ILC7 and 
NILC7 map obtained by the Minkowski functionals for the 
surrogates Arrange [120,300], which is covering the first 
peak of the power spectrum. The three Minkowski function- 
als show indications for phase correlations for both maps. 
Yet the directions of these deviations on the sky totally dif- 
fer from each other. Since the two compared maps differ in 



their resolution we extend the tests with a comparison of 
the NILC7 at a resolution decreased to one degree FWHM 
(not shown here). We generate the corresponding surrogate 
maps of first and second order and calculate the statistics 
for the Minkowski functionals and the scaling indices. The 
signatures for phase correlations detected by the Minkowksi 
functionals differ from the results of the fully resolved NILC7 
map and they do not resemble the signatures of the ILC7 
map. For the SIM there is more concordance between the 
findings of the differing resolutions of NILC7. We claim that 
the reason for disparities between the ILC and NILC maps 
on intermediate length-scales, 120 ^ / ^ 300, is not only the 
difference in beam resolutions but even more the difference 



in the foreground cleaning. As explained inlRath et al. ( 2011 ) 



foreground removal can induce phase correlations, especially 
on these intermediate ^-ranges. The detected phase correla- 
tions in the A^-range of [120, 300] depend on the foreground 
cleaning, the map resolutions and the response of the image 
analysis techniques. If there are additional deviations from 
Gaussianity with cosmological origin in the maps on these 
^-ranges we cannot distinguish them from systematics so far. 
Our findings clearly show that the found non-Gaussianities 
in WMAP data are scale-dependent and can have different 
origin. 

So far, we have analysed the deviations from Gaussian- 
ity in areas of apex angle ix around the sky which leads to 
a solid angle of 27T, so a full hemispherical analysis. We de- 
crease the apex angle to ty/2 and calculate the Minkowski 
functionals and scaling indices for areas with a size of only 
30% (solid angle 0.67r) of the former hemispheres. The cor- 
responding 5-maps are shown in Figure [8] for the x 2 statistic 
of the three Minkowski functionals and the SIM. The tt/2 
study analyses smaller sky patches of the CMB. It can fea- 
ture a better spatial localisation of the phase correlations. 
However, a smaller number of pixels might constrict the 
detection of a weaker signal. Comparing between the two 
classes of surrogates on smaller areas of the sky maps, we 
find no signal for no n- Gaussianity in the northern ecliptic 
sky, whereas we detect individual spots in the southern sky 
which indicate phase correlations in these parts of the sky. 
The fact that we do not detect these spots in the northern 
sky makes ecliptic systematics from observations less likely. 
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olution parameter of the maps is decreased to N s id e = 64 
which does not change the results for low-^ ranges. One can 
immediately see that the ^-values of the /NL al_s i mma tions 
are nearly always smaller compared to the original data set. 
The few exceptions where the signal of the simulations lies 
close to the NILC7 data feature a value of f^£ al = ±1000. 
These extreme values are already ruled out by recent analy- 
ses that resulted in current constrai nts for | f^ cal I well belo w 
100 (e.g. /^°£ al = 32 ±21 (68% CL), |Komatsu et al.| ( [20lT] )) . 
For the scaling indices, we obtain similar results, which are 
shown in Figure [l0| For the maximum plot, only one simu- 
lation lies above the results for the NILC7 map, which is at 
/^ al = 100. For the mean, we obtain three points with a 
slightly higher value than the data, and none for the stan- 
dard deviation. 

One has to conclude that the CMB simulations with 
/NL al cannot reproduce the low-^ anomalies we found in the 
WMAP data. This means that the detected large-scale NGs 
and asymmetries in the data do not correspond to the type 
of NG which is described by /nl^ but stem from a differ- 
ent origin. The search for the source of the deviations from 
Gaussianity has to be continued in future analyses. 
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Figure 10. Maximum, mean and and standard deviation of the 
absolute values of S(x 2 )a of ten f^ al simulations compared to 
one NILC7 original map (red line). The HEALPix resolution of 
the maps is N s ^ e = 64. The ^-range for the method of the surro- 
gates is A£ = [2,20]. 



However, it cannot be ruled out that the detected spots are 
correlated with unknown foregrounds. Further study on the 
origin of the found spots and comparison with other detected 
anomalous spots in the CMB is needed. 



5.2 /NL-dependent Simulations 

The analysis of the ILC and NILC7 surrogate maps with 
the Minkowski functionals as well as the scaling indices re- 
veals HOCs in the original maps. In order to test whether 
these phase correlations can be reproduced by /NL al_m °dels 
we analyse constrained realisations of the temperature 
maps with varying /^l^ -parameter with values /nl^ 1 = 
[0, ±100, ±1000] by means of the surrogate method for A£ = 
[2,20]. Figure [9] shows the results for the corresponding 
S(x 2 ) statistics for the three Minkowski functionals. We plot 
the maximum, mean and standard deviation of S(x 2 ) f° r ten 
simulations and the NILC7 original map. The HEALPix res- 



6 CONCLUSIONS 

In addition to our previous work based on the use of sur- 
rogate maps we analysed latest WMAP experiment CMB 
maps with respect to asymmetries and scale-dependent 
non- Gaussianity. The surrogates are generated by a scale- 
dependent shuffling of Fourier phases while all other prop- 
erties of the maps are preserved. In this work we focus on the 
Minkowski functionals calculated for the ILC7 and NILC7 
maps as a scale-independent measure being sensitive on the 
HOCs of the maps. We compare these new results with 
the scale-dependent scaling index method calculated for the 
NILC7 map in this work and for the ILC7 map in previ- 
ous works. We find that both measures detect highly signif- 
icant signatures for phase-correlations and therefore devia- 
tions from Gaussianity, and furthermore ecliptic hemispheri- 
cal asymmetries for the interval A£ — [2, 20] in both ILC and 
NILC maps. The reduction of the analysed sky region down 
to 30% of the former hemispheres shows no signal for NG in 
the northern ecliptic sky. In the south we detect individual 
spots of NG. 

If the findings are indeed of intrinsic nature they would 
disagree with the predictions of isotropic cosmologies with 
single-field slow-roll inflation. 

The two different image analysis techniques reveal very 
consistent results on the \ow-£ range for both maps. The 
signatures we find for A£ = [120, 300] show less agreement 
between ILC and NILC and differ from the low-^ range re- 
sults. They also depend on the image analysis technique and 
the resolution of the map. As discussed in our earlier works 
NGs on these scales can be induced by foreground cleaning 
and can be subject to secondary anisotropics. 

The constrained realisations of the CMB with varying 
/NL al cannot parametrise the detected NGs and asymmetries 
on largest scales. /nl &1 might still be a suitable parametrisa- 
tion for smaller scales (larger-^ ranges). Also, other types of 
/nl, as the equilateral and orthogonal type, also the #nl pa- 
rameter could characterise the found NGs. A deeper study 
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Figure 9. Maximum, mean and standard deviation of the absolute values of S(x 2 )m* of ten f^ al simulations compared to the NILC7 
original map (red line) for area (left column), perimeter (middle) and Euler (right). The HEALPix resolution of the maps is N s [d e = 64. 
The ^-range for the method of the surrogates is A£ = [2, 20]. 



of different inflationary models, as for example Bianchi type 
models, that represent appropriate test candidates for the 
found anomalies, is required. 

Future investigations, e.g. of upcoming data of the 
Planck satellite, will shed more light on the open questions 
regarding instrumental constraints, observational systemat- 
ics, map making influences and resolution problems. 
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